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ABSTRACT: 

\n  this  paper  an  efficient  finite  element  method,  which  accounts  for  the 'translating 
singularities  near  the  tip  of  a  dynamically  propagating  crack,  is  presented.  Two 
aspects  of  analyses  of  dynamic  fracture,  namely,  (i) determination  of  dynamic 
stress-intensity  factors  from  given  crack  velocity  history,  and  (ii)determi- 
nation  of  crack-velocity  history  (and  arrest)  from  given  dynamic-f racture-tough- 
uess  versus  crack-velocity  data  are  considered.  Results  are  presented  and  dis¬ 
cussed  for  wedge-loaded  rectangular,  as  well  as  tapered,  double-centilever  beam 
specimens.  .. 
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Dynamic  Fracture,  Propagating  Crack  Eigen  Functions,  Translating  Singularities, 
Finite  Elements,  Dynamic  Stress-intensity  Factor,  Dynamic  Fracture  Toughness, 
Prediction  of  Crack  Propagation. 

INTRODUCTION 

A  number  of  analytical  and  experimental  investigations,  into  dynamic  propagation 
and  arrest  of  cracks  in  structural  metals,  have  been  conducted  in  recent  years 
due,  in  part,  to  the  important,  technical  problems  in  asessing  the  structural  in¬ 
tegrity  of  flawed  ship  structures  and  of  nuclear  pressure  vessels,  and  the  need 
to  develop  criteria  for  crack  arrest  in  these  situations.  Analytical  solutions 
exist,  for  the  most  part,  for  boundary-value  problems  in  linear  elastodynamics  for 
cracks  in  unbounded  media.  The  problem  of  finite  bodies,  wherein  the  interaction 
of  stress-waves  reflected  from  the  boundaries  and/or  the  other  crack-tip  with  the 
considered  crack-tip,  play  an  important  role,  is  amenable  for  analysis,  in  most 
cases,  only  though  numerical  techniques.  The  literature  on  numerical  methods, 
such  as  finite  elements  and  finite  differences,  in  dynamic  fracture  has  been 
succinctly  reviewed  by  Kanninen  (1978).  Most  of  the  techniques  reviewed  by 
Kanninen  (1978)  fall  into  the  category  of  "node-release  techniques". 

Recently,  Atluri,  Nishioka,  ar  '  Nakagaki  (1979),  and  Nishioka,  and  Atluri  (1980 
a,b,c)  have  presented  formulations  for  a  moving-singularity  finite  element  method 
for  simulating  fast  crack  propagation  in  finite  bodies.  An  energy-consistent 
variational  statement  is  first  developed  (Atluri,  Nishioka,  and  Nakagaki,  1979; 
Nishioka,  and  Atluri,  1980  a),  as  a  basis  for  the  above  "moving  singular-element" 
procedure  of  dynamic  fracture  analysis.  The  above  numerical  procedure  has  been 
applied  (Nishioka  and  Atluri,  1980  b,c)  in  studying  the  finite-domain  counterparts 


of  the  problems  with  unbounded  domains  for  which  analytical  solutions  exist  (Bro- 
berg,  1960;  Baker,  1962;  Nilsson,  1972;  Freund,  1973).  The  obtained  numerical  re¬ 
sults  have  been  found  to  correlate  excellently  with  the  cited  analytical  solutions 
for  infinite  domains,  during  the  time  intervals  for  which  these  analytical 
solutions  may  be  considered  as  valid.  The  computed  solutions  beyond  these  times, 
and  the  knowledge  of  the  various  time  periods  involved  for  wave-interaction  in 
finite  bodies,  indicated  both  qualitatively  as  well  as  quantitatively  the  effects 
of  stress-wave  interactions  on  dynamic  K-factors  for  cracks  propagating  in  finite 
bodies. 

In  this  paper,  attention  is  primarily  focused  on  the  "application  phase",  viz., 
the  determination  of  crack-tip  velocity  history  (and  possible  crack  arrest)  from 
a  given  dynamic- fracture- toughness  (K^)  versus  the  crack-velocity  (Z)  relation, 
using  the  above  described  numerical  procedure.  To  this  end,  the  experimental 
results  of  Kalthoff,  Beinert,  and  Winkler  (1977,  1978)  on  wedge-loaded  rectangu¬ 
lar  (WL-RDCB)  as  well  as  tapered  (WL-TDCB)  double  centilever  beam  specimens  are 
considered.  In  the  first-phase  of  the  presently  reported  study,  the  experimental 
results  for  Z  versus  time  for  the  WL-RDCB  specimen  are  simulated  so  as  to  numeri¬ 
cally  determine  the  dynamic  K-factor  variation  with  time  and  compare  it  with  the 
reported  experimental  (caustics)  measurements  of  dynamic  K-factors.  In  the  second 
phase  of  the  study,  using  the  reported  (Kalthoff,  Beinert,  and  Winkler,  1977;  1978) 
experimental  results  for  K^D  versus  Z  relations,  and  the  crack- initiation  stress- 
intensity  factor  K  ,  the  crack-propagation  history  in  WL-RDCB,  and  WL-TDCB  speci¬ 
mens  is  numerically^predicted  and  compared  with  cited  experimental  results.  The 
oresent  results  are  also  compared,  where  possible,  with  independent  numerical 
results  of  Kobayashi  (1979).  Discussions  are  presented  on  the  possible  specimen- 
geometry  dependence  of  the  K.^  vs  Z  relation. 

OUTLINE  OF  TRAN SLAT ING- S INGULARITY  FINITE-ELEMENT  METHOD 

As  derived  by  Atluri,  Nishioka,  and  Nakagaki  (1979),  the  virtual  work  equation 
applicable  to  a  system  of  finite  elements  used  to  model  the  cracked  body,  be¬ 
tween  the  times  t^  and  t^,  during  which  a  Mode-I  type  crack-growth  by  an  amount 
AZ  takes  place  can  be  written  as: 


In  the  above,  V  is  the  domain  of  the  singularity-element;  Z  is  the  crack-sur¬ 
face  contained  within  the  boundary  of  V  at  time  t.;  the  suplr-scirpt  (+)  refers 
to  the  upper-half  of  the  crack-surface  which  only  is  modelled  in  symmetric 
M?fie-I  case;  AZ  is  the  amount  of  crack  growth  between  times  t^  and  t  ;  V  is  the 
n  regular  (isoparametric)  finite  element  in  the  cracked-body;  S  is  that  por¬ 
tion  of  the  boundary  (away  from  the  carck)  of  ;  o,  T,  u,  <$e,  5u  are,  respecti¬ 
vely,  the  stress  tensor,  applied  traction  vector,  displacement  vector,  virtual 
stain  tensor,  and  virtual  displacement  vector;  the  superscripts  (1)  and  (2) 
denote,  respectively  the  quantities  at  times  t.  and  yl  is  the  unit  normal 
vector  to  AZ.  In  the  above  the  notation  g: 6e=ojj 5e^j  and  T.6u=T^6u^  has  been  used. 
Note  that  at  time  t»,  when  the  crack  has  opened  Dy+an  additional  amount  AZ,  the 
admissible  virtual  displacements  are  such  taht  <5uj  1*0  at  AZ. 
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In  the  present  procedure,  the  displacement  field  within  the  singular-element,  V  , 

s 

is  assumed  to  be: 

u  (£,x„,t)  =  Eu*  (£,x  v)B  (t)  +  Rigid  body  modes  (a=l,2)  (2) 

a  2  m  am  2  m 

When  u  are  eigen-function  solutions  for  the  elastodynamic  problem  of  steady- 
state  ^with  respect  to  an  observer  moving  with  the  crack-tip)  crack-propagation 
at  constant  velocity  v,  as  given  in  Atluri,  Nishioka,  and  Nakagaki  (1979).  If  x^ 
(along  the  crack-axis) and  x^  (normal  to  crack-axis)  are  fixed  cartesian  coordi¬ 
nates,  (£,x  )  is  the  coordinate  system  that  moves  with  the  crack-tip,  ie., 

£=x^-vt.  In  the  present  procedure,  the  assumed  displacement  field  at  times  t^ 
andt_,  respectively,  is  taken  to  be  that  corresponding  to  Eq.  (2)  wherein  the 
velocity  v  takes  on  the  values  v^  and  v^,  viz.,  the  velocities  at  times  t^  and  t 
It  is  further  seen  from  the  properties  of  the  above  eigen-functions  (Atluri, 
Nishioka,  and  Nakagaki,  1979)  that  8^(t)  in  Eq.  (2)  is  in  fact  the  dynamic  stress- 
intensity  factor  and  hence  may  be  computed  directly. 


The  total  velocity  and  acceleration  of  a  material  particle  in  V  are  then: 

it  *  it  m  ® 

u  =  Eu  8  -  vEu  g 

a  m  am  m  m  am,£  m 

if  ••  it  ,  0  it 

and  li  =  Eu  8  -  2vEu  r8  +  v  Eu  8 

a  m  am  m  m  am,  t,  m  m  am,t,t,  m 


(3) 

(4) 


where  (  ),r=3(  )/3£.  The  compatability  of  displacement,  velocity,  and  acceleration 
at  the  boundary  of  V  with  the  surrounding  regular  elements  is  satisfied  through 
a  least-squares  technique  (Atluri,  Nishioka,  and  Nakagaki,  1979). 


Assuming  that,  between  times  t  and  t^>  crack  growth  occurs  (which  can  be  deter¬ 
mined,  as  explained  below,  in  the  so-called  "application  phase"  calculations 
using  the  given  material  characteristic  of  a  K  _  versus  E  relation;  or  is  known 
apriori,  in  the  so-called  "generation  phase"  calculation,  ie. ,  in  the  case  of 
simulation  of  given  E  versus  time  data),  the  singularity  element  is  translated,  by 
an  appropriate  distance  AE  from  its  location  at  time  t^  as  shown  in  Fig.  1. 


FIGURE  1 
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Note  that,  unlike  in  the  node-release  technique  of  Kobayashi  (1979),  the  above 
shown  amount  AE  is  not,  in  any  way,  related  to  the  distance  between  any  two  ad¬ 
jacent  finite  element  nodes  at  time  t.^.  Note  also  from  Fig.  1  that,  as  the  singu¬ 
lar-element  is  translated,  the  surrounding  regular  elements  are  appropriately 
distorted.  After  a  sufficient  amount  of  translation,  say  M.(AE),  of  the  singular- 
element,  the  distortion  of  the  regular-elements  may  become  accute.  In  such  a  case, 
the  regular  finite  element  mesh  is  readjusted,  as  described  in  Nishioka,  and 
Atluri  (1980a). 


The  finite  element  equation  obtained  by  using  the  variational  principle  of  Eq. 


(1)  can  be  shoen  to  be: 

m 2  +  ~s2  =  92+9x  -  for  VR 

?s-s2  +  ?sSs2  +  =sqs2  =  9S  for  Vs 


(4) 


where  V  is  the  region  spanned  by  regular  elements.  It  is  shown  by  Nishioka  and 
Atluri  F1980a)  that  out  of  all  the  matrices  appearing  in  Eq.  (4),  two,  viz.,  Ks 
and  D  are  unsymmetric,  due  to  the  use  of  the  eigen-functions  in  the  moving 
coordinate  system  (S.x.)  as  in  V  .  The  Equs.  (4)  are  numerically  sovled  in  the 
present,  using  the  well-known  Newmark-B  method.  Once  the  displacements,  veloci¬ 
ties,  and  accelerations  (viz.,  9s2»9s2»  an<*  9s2)  t^e  nodes  of  V  are  solved 
for  from  Eq.  (4),  the  undetermined  parameters  g  ,  §  ,  and  B  ,  (andshence  the  dy¬ 
namic  stress-intensity  factor,  and  its  first  anl?  second  time-derivatives)  can  be 
computed  directly. 


CALCULATION  OF  K  (t)  FROM  GIVEN  Z(t)  DATA  ("GENERATION  PHASE"): 

Since  the  crack-velocity  v=£(t),  is  known,  using  the  velocitites  v^  and  v^  at  times 
t^  and  t^,  respectively,  the  finite  element  equations  (4)  can  be  easily  construc¬ 
ted.  From  these,  as  explained  above,  B^  (which  is  in  fact  equal  to  K  ),  B^,  and 
B^  can  be  directly  computed  at  time  t 2  as  well  as  at  all  subsequent  times. 

PREDICTION  OF  Z(t)  FROM  GIVEN  VERSUS  E(t)  RELATION  ("APPLICATION  PHASE"): 

Let  the  solution  upto  time  t^  be  assumed  to  be  known.  In  order  to  find  the  solu¬ 
tion  at  t2  from  Eq.  (4),  the  crack-velocity  at  t2,  ie. ,  v=E(t  )  must  be  found. 

To  this  end,  it  is  first  noted  that  the  dynamic  stress-intensity  factor  can  be 
written  as: 


K  =  K-^t.v). 


(5) 


Since,  in  the  present  procedure,  the  velocity  of  crack-propagation  is  assumed  to 
be  constant  within  each  time-step,  an  approximate  procedure  to  predict  the  ve¬ 
locity  at  [tj+(At/2)]  will  be  sought.  Using  double  Taylor  series  expansion,  it 
is  seen  from  Eq.  (5)  that: 


KiP(ti+  r;VAv)  =  I=o  bJ(r)  h +  Av  hl\(t i*vi) 


(6) 


where  KT  is  the  predicted  K  at  t1+(At/2)  while  K_(t-,v-)  is  known.  One  can  now 
rewrite  pEqn.  (6)  as: 

2  . 

KIp  =  KI(t1,v1)  +  (fL)KI(t1,v1)  +  KI(t1,v1)  +  R 

•  2** 

a  B1(t1)  +  (|t)B1(t1)  +  lAll  B1(t1)  +  Rs  K*p  +  R  (7a, b) 

where,  (  )=3(  )/3t,  R  is  the  "residue"  (all  the  omitted  terms  from  the  right- 
hand  side  of  Eq.  (6)),  and  as  explained  above,  in  the  present  procedure, 

Kl(tl)=3l(Cl>,  etc* 


Using  the  predicted  K-  of  Eq.  (7),  and  the  given  K  versus  £(t)  relation,  the 
crack-velocity  v(or  £/Pat  time  [t  +(At/2)]  can  be  predicted.  Using  this  crack- 
velocity  value,  the  finite  element  system  of  equation  (4)  at  time  t_  can  be 
constructed,  and  from  these,  the  actual  dynamic  stress-intenstiy  factor  K (t?) 
can  be  computed.  Thus,  the  actual  at  t^+(At/2)  can  be  computed,  as: 

KI(tl+2t)  +Kj(t2)] 


(8) 


The  correlation  between  the  predicted  K  of  Eq.  (7)  and  the  actual  K  of  Eq.  (8) 
can  be  seen  to  depend  on  the  residue  'R+pof  Eq.  (7).  To  ensure  this  correlation, 
a  further  approximation  is  introduced  in  the  present  work  that  the  residue  R  at 
time  tj+(At/2)  can  be  approximated  by  its  known  value  at  t..-(At/2),  in  the  generic 
case.  Thus,  in  the  present  procedure,  Eq.  (7)  is  re-approximated  as: 

Vc+  r>  -  W  +  +  -  tVv  f>-  Arici- 

In  the  present  calculation,  when  Eqn.  (7)  with  R=0  was  used  an  error  of  the 
order  of  3 %  between  Kj  and  Kj  was  noted.  However,  when  Eq.  (9)  was  used,  this 
error  reduced  to  the  o?der  of  0.5%.  Further  algorithmic  details  of  the  predic¬ 
tion  phase  along  with  detailed  error  analysis  studies,  which  are  omitted  here, 
will  be  reported  elsewhere. 

RESULTS  OF  GENERATION  CALCULATIONS: 

The  presently  used  finite  element  mesh  for  the  upper  half  of  WL-RDCB  specimen 
(Kalthof f  and  collegues,  1977)  is  shown  in  Fig.  2.  Note  that,  in  the  present 
translating  singular-element  procedure,  this  mesh  represents  the  configuration 
at  the  beginning  of  the  crack-propagation.  The  singular-element  is  shown  by  the 
shaded  area  in  Fig.  2.  Plane-strain  conditions  are  invoked.  The  experimentally 
measured  crack-length  versus  time  history  (as  well  as  the  crack-velocity  versus 
time  history)  for  the  WL-RDCB  specimen  (reported  by  Kalthof f,  Beinert,  and  Winkler, 
1977)  with  a  crack-initiation  stress-intensity  factor  K  (larger  than  fracture 
toughness  K  )of  2.32  MN/m3/2  are  shonw  in  Fig.  3.  q 
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FIGURE  2 


FIGURE  3 


Also  shown  in  Fig.  3  are:  (£)  the  presently  calculated  dynamic  K-factor,  based 
on  (a)  direct  evaluation  of  as  explained  earlier,  (b)  using  the  crack-opening 
energy  as  computed  from  the  present  variational  formulation,  and  the  crack-tip 
integral,  and  (c)  using  the  crack-opening  energy  as  computed  from  a  global 


v  [m/sec] 


energy  balance,  and  ( ii )  the  experimentally  determined  dynamic  K-factors  and 
{Hi)  the  numerical  results  of  Kobayashi  [1979].  It  is  seen  that  the  present 
results  under  (i) a  and  b  are  almost  identical.  It  is  also  seen  that  the 
present  results  correlate  excellently  with  those  of  Kalthoff  and  co-workers 
(1977)  during  the  first  (2/3)  of  the  propagation  history.  However,  during  the 
last  of  (1/3)  of  propagation  history,  the  present  results  exhibit  a  pronounced 
peak  as  compared  to  the  experimental  results.  It  is  felt  that  this  may  be  partly 
due  to  the  fact  that  the  present  is  a  linear  elastodynamic  analysis  while  the  ex¬ 
perimental  material  (Araldite  B)  is  viscoelastic  in  nature.  It  is  interesting  to 
note  that  the  numerical  results  of  Kobayashi  (1979)  also  exhibit  a  pronounced 
peak  as  in  the  present  results,  except  that  this  occurs  much  earlier  in  Kobayashi' s 
results.  The  problem  was  reanalysed  by  slightly  perturbing  the  data  for  E(t) 
curve.  These  results  (not  presented  here  for  want  of  space)  indicate  that  these 
slight  perturbations,  while  not  affecting  the  K-f actor  solution  during  the  first 
(2/3)  of  propagation,  do  alter  the  magnitude  of  the  above  mentioned  peak  in  the 
latter  (1/3)  of  propagation.  These  results  indicate  the  need  for  the  extreme 
care  with  which  experimental  data  for  Z(t)  should  be  recorded. 


Results  of  Application  Phase: 


The  experimental  (Kalthoff  and  co-workers,  1977;  1978)  results  for  K^  versus 
L  relation,  for  WL-RDCB  and  WL-TDCB  specimens,  that  are  used  in  the  present  anal¬ 
yses  are  shown  in  Fig.  4. 


3/2 

Two  RDCB  specimens,  with  K^  values  of  2.32  and  1.33  MN/m  ,  respectively,  were 
analyzed  using  a  plane-strain  assumption.  Using  the  KT_  versus  £  relation  of  Fig. 
4,  the  propagation  histories  in  both  specimens  were  predicted,  using  the  earlier 
described  procedure.  The  results  for  the  two  RDCB  specimens  are  shown  in  Figs. 5 
and  6  representively ,  along  with  comparison  experimental  results  of  Kalthoff  and 
co-workers  (1977)  and  available  numerical  results  of  Kobayashi  (1979).  It  is 
seen  that  the  present  predictions  agree  excellently  with  the  cited  experimental 


\ 


results.  However,  it  is  seen  that,  for  the  RDCB  specimen  with  =1. 33MN/m^/2 > 
the  predicted  crack  arrest  length  is  about  19%  higher  than  the  experimental 
lenghth.  Once  again,  it  is  believed  that  this  descrepancy  is  due  to  the  inade¬ 
quacy  of  the  present  linear-elastic  material  modeling.  The  variation  of  dif¬ 
ferent  energy  quantities  of  interest  during  this  predicted  crack-propagation 
history,  in  the  first  RDCB  specimen  (K  =i. 33MN/m^' ^)  is  shown  in  Fig.  7.  Since 
each  of  the  indicated  quantities  T,U,  and  F,  were  calculated  independently 
using  the  present  procedures,  the  fact  that  their  sum  correlates  excellently  with 
the  input  energy  W  lends  credence  to  the  present  calculations. 
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FIGURE  7 


The  finite  element  model  for  TDCB  specimen  (at  initiation  of  propagation)  is 
shown  in  Fig.  8.  A  value  of  2.08MN/m^^  is  used.  The  predicted  crack-propa¬ 
gation  history,  based  3n  the  K  vs  E  relation  of  Fig.  4,  is  shown  in  Fig.  9. 

Once  again  good  correlation  of  the  present  results  with  the  experimental  ones 
(Kalthoff  and  co-workers,  1978)  is  noted. 


Results  were  also  obtained  in  a  study  of  the  possible  specimen  dependency  of  the 
versus  E  relation.  These  results,  along  with  details  of  prediction  phase  cal¬ 
culations  for  several  other  specimens,  which  are  omitted  here  for  space  reasons, 
are  mentioned  in  the  oral  presentation  of  this  paper  and  will  be  documented  else¬ 
where. 
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